Proc. 29th International Cosmic Ray Conference, Pune, India (2005) 



Thermal Particle Injection in Nonlinear Diffusive Shock Acceleration 

Donald C. Ellison'^, Pasquale Blasi'', Stefano Gabici'^ 

(a) Physics Dept., North Carolina State Univ., Raleigh, NC 27695, U.S.A. 

(b) High Energy Astrophysics Group, INAF, Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, 1-50125, Firenze, Italy 

(c) Humboldt Fellow, Max-Planck-lnstitut fuer Kernphysik, Saupfercheckweg 1, 691 17 Heidelberg, Germany 
Presenter: Don Ellison (don_ellison@ncsu.edu), usa-ellison-D-abs3-ogl4-oral 



Particle acceleration in collisionless astrophysical shocks, i.e., diffusive shock acceleration (DSA), is the most 
likely mechanism for producing cosmic rays, at least below 10^^ eV. Despite the success of this theory, several 
key elements, including the injection of thermal particles, remains poorly understood. We investigate injec- 
tion in strongly nonlinear shocks by comparing a semi-analytic model of DSA with a Monte Carlo model. 
These two models treat injection quite differently and we show, for a particular set of parameters, how these 
differences influence the overall acceleration efficiency and the shape of the broad-band distribution function. 

1. Introduction 

We compare a recent semi-analytic (S A) model of non-linear diffusive shock acceleration (DSA) f2l f3l f4l with 
a well-established Monte Carlo (MC) model (e.g., 1 12|| 15|| 11 1). Both include a thermal leakage model for 
injection and the non-linear (NL) backreaction of accelerated particles on the shock structure, but they do these 
in very different ways. Also different is the way in which the particle diffusion in the background magnetic 
turbulence is modeled. Our limited comparison shows that the important NL effects of compression ratios ^ 4 
and concave spectra do not depend strongly on the injection model as long as injection is efficient. A fuller 
understanding of the complex plasma processes involved, particularly if injection is weak, will require particle- 
in-cell (PIC) simulations (e.g., 1 13 1), but these simulations, which must be done fully in three-dimensions 1 14|, 
cannot yet be run long enough, in large enough simulation spaces, to accelerate particles from thermal to 
relativistic energies in order to show strong NL effects. For now, approximate methods must be used. 

The Monte Carlo model is more general than the semi-analytic model,' but it is considerably slower compu- 
tationally. Since it is important in many applications, such as hydro models of supernova remnants, to include 
the dynamic effects of nonlinear DSA in simulations which perform the calculation many times (e.g., |8|), a 
rapid, approximate calculation, such as the SA one discussed here, is useful. 

2. Models 

The main features of the two models that are relevant for our comparison here are: 

Injection: In the semi-analytic model, a free injection parameter, ^sA, determines the fraction of total particles 
injected into the acceleration mechanism and the injection momentum, pmj ■ Specifically, pmj = CsA Pth, where 
Pth — \/2mpkTY)s (Tds is the downstream temperature). The fraction, rjsA, of unshocked particles crossing 
the shock which become superthermal in the SA model is rjsA — (4/37r'^/^)(rsub — 1)Csa^~^^*' where rgub 
is the subshock compression ratio. The fraction rjsA, which is approximately the number of particles in the 
Maxwellian defined by Tds with momentum p > pinj, is determined by requiring the continuity, atpinj, of the 
Maxwellian and the superthermal distribution (see H). Since pth depends on the injected fraction, the solution 
must be obtained by iteration. 

In the Monte Carlo model, the injection depends on the scattering assumptions. We assume that particles 
pitch-angle scatter elastically and isotopically in the local plasma frame and that the mean free path (mfp) 



^For instance, the MC model can treat a specific momentum dependence for the scattering mean free path, pailicle acceleration in 
relativistic shocks l9l 1 101 . and NL effects in oblique shocks ITTI . 
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is proportional to the gyroradius, i.e., A oc r^, where = pc/{qB). With these assumptions, the injection 
is purely statistical with those "thermal" particles which manage to diffuse back upstream gaining additional 
energy and becoming superthermal. Note that in this scheme the viscous subshock is assumed to be transparent 
to all particles, even thermal ones, and that any downstream particle with v > U2 has a chance to be injected. 
Here, U2 is the downstream flow speed. 

For comparison with the SA model, we have included an additional parameter, ii^J^g, to limit injection in the 
MC simulation. Only downstream particles with v > v^^^^ are injected, i.e., allowed to re-cross the shock into 
the upstream region and become superthermal. In our previous Monte Carlo results, with the sole exception of 
we have taken I'J^^g = 0. 

Momentum dependent diffusion: Diffusion is treated very differently in the two models. As just mentioned, 
the MC simulation models pitch-angle diffusion by assuming a X{p) (see 1 10 1). The semi-analytic model does 
not explicitly describe diffusion but assumes only that the diffusion is a strongly increasing function of particle 
momentum p so that particles of different p interact with different spatial regions of the upstream precursor 
Eichler |5 1 used a similar procedure. With this assumption, particles of momentum p can be assumed to feel 
some average precursor fluid speed Up and an average compression ratio ~ Up/u2 (see |4J). 



Uo=10'* km/s, p„„^=100 rripC, Ms=200, M^=229 
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Figure 1. Phase space distributions f{p) times for a MC model (solid curves) and a SA model (dashed curves). In the 
left panel w^J^^ = 0, and in the right panel — 4it2 in the MC model. Here, U2 = uo/r is the downstream plasma 

speed in the shock frame, Mo is the shock speed, and the Mach numbers, compression ratios, and shocked temperatures are 
indicated. In all cases, the shock is parallel and Alfven heating is assumed in the shock precursor (TJ. 

The different way diffusion is treated influences not only injection, but also the shape of the distribution func- 
tion f{p) [f{p) is the momentum phase space density, i.e., particles/(cm'^-d"^p)]. Both models give the charac- 
teristic concave f{p) which hardens with increasing p and, since the overall shock compression ratio r can be 
greater than 4, this spectrum will be harder than p^* at ultra-relativistic energies. In the results we show here, 
the acceleration is limited with a cutoff momentum so / (p) cuts off abruptly at Pmax- More realistic models will 
show the effects of escape from some spatial boundary (e.g., finite shock size) or from a finite acceleration time. 
In either case, the spectrum will show a quasi-exponential turnover, e.g., f{p) oc p~'^ exp [— a"^ (p/pmax)"] 
f7). where a is included to emphasize that the detailed shape of the turnover depends on the momentum de- 
pendence of the diffusion coefficient near Pmax- 
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Thermalization: There is no thermalization process in the MC simulation in the sense of particles exchanging 
energy between one another because particles scatter elastically in the local frame. However, a quasi-thermal 
low energy distribution is created as unshocked particles cross the shock from upstream to downstream at 
different angles and receive different fractions of the speed difference, uq — U2- For the parameters used here, 
the low-energy peak is essentially a Maxwellian when wj^^g — 0. 

In the semi-analytic model, the shocked thermal pressure and density are determined from the conservation 
relations and these are translated to a Maxwell-Boltzmann distribution. 
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Figure 2. Low energy portions of the spectra shown in Fig. Q where the solid histograms are the MC results and the 
dashed curves are the SA results. The heavy dotted histograms show the distribution of particles that were injected in the 
MC model. The solid vertical lines are drawn at pinj, the momentum at which particles are injected in the SA model. For 
t'ftSs = 0, r?MC = 5.7 X 10"^ and77sA = 4.6 x 10"^ while for Wt^J^^, = 4^2, r^MC = 4.1 x 10"^ andr;sA = 2.5 x 10~^ 
The light-weight dotted curve shows the Maxwellian (Tds = 2.2 x 10^ K) that would have resulted if no DSA occurred. 

3. Results 

In Fig.^we show distribution functions (times p^) for a set of parameters and with = (left panel) and 

^thros ~ '^^s (right panel). In both panels, the solid curves are the MC results and the dashed curves are the 
SA results where ^sA has been chosen to provide the best match to the MC spectra.^ In Fig.|2lwe show the 
low energy portions of the spectra along with the distributions of injected particles in the MC model (dotted 
histograms). The vertical solid lines are at pinj and they show the transition between the Maxwellian and the 
superthermal population in the SA results. 

The important aspects of the plots are: 

(i) The broad-band match between the two very different calculations is quite good, particularly for w^J^g — 0. 
Both models show the important characteristics of nonlinear DSA, i.e., r ^ 4, concave spectra, and a sharp 
reduction in the shocked temperature from test-particle values. 

(ii) The shocked temperature depends on v^^^ with weaker injection (i.e., larger f JflJ^g) giving a larger Tds.'' 



■^Note that once the input parameters (i.e., mq, Pmax, and the sonic and Alfven Mach numbers: A/g and Ma) are fixed, the 

only parameter varied to match the spectra shown in Figs.0and|2]is 5sA- The normalization is not adjusted. 

•^The light-weight dotted curve in the left panel of Fig.|2|is a Maxwellian at Tds = 2.2 x 10^ K, the temperature of the shocked gas 
without DSA. 
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(iii) The distribution f{p) is harder near Pmax in the MC results than with the SA calculations. A common 
prediction of SA models is that the shape of the particle spectra at Pmax has the form p^3.5 ^j^g shock is 
strongly modified and the diffusion coefficient grows fast enough in momentum. This can be demonstrated by 
solving the equations in the extreme case of a maximally modified shock and approximating the spectrum with 
a power law at high momenta. The MC model makes no such power-law assumption. 

(iv) The minimum in the p^f{p) plot occurs at p > mpC in both models. The transition between f{p) softer 
than p^'^ and f{p) harder than varies with shock parameters and increases as Pmax increases. This is an 
important difference from the algebraic model of Berezhko & Ellison[JJ where the minimum is fixed at m^c. 

(v) For the particular parameters used in these examples, the overall compression ratio is relatively insensitive 
to I'JjJ^g, but shocks having other parameters may show a greater sensitivity. Note that Alfven wave heating 
is assumed in all of the results presented here. If only adiabatic heating was assumed, the compression ratios 
would be much higher (see (TJ). 

(vi) In all cases, the MC model injects more particles than the SA model but the average energy of the injected 
particles is less, as indicated by the peak of the curve labelled 'MC inj' vs. the 'SA inj' energy in Fig.|2] 

(vii) In contrast to v^^^ = 0, the 'thermal' part of /(p) with v^^^ = 4u2 (Fig.|2j shows large differences 
in the two models. While both conserve particle, momentum, and energy fluxes so that the broad-band f{p) 
matches well for a wide range of parameters, the different treatments of the subshock lead to large differences 
in the critical energy range 2 < [i?/(/cTDs)] ^ 5. This offers a way to distinguish these models observationally. 

4. Discussion and Conclusions 

Using two approximate acceleration models, we have shown that the most important features of NL DSA, i.e., 
r ^ 4 and concave spectra, are robust and do not strongly depend on the injection model as long as injection 
is efficient. If injection is weak, as might be the case in highly oblique shocks, accelerated spectra will depend 
more on the details of injection, at least in the transition range between thermal and superthermal energies. 
Also, the relative efficiencies for injecting and accelerating electrons vs. protons or protons vs. heavier ions 
may require a more detailed description of injection, as may be provided by future PIC simulations. 

Acknowledgements: D.C.E. wishes to acknowledge support from a NASA grant (ATP02-0042-0006) and 

5. G. acknowledges support from the Humboldt foundation. 

References 

[1] Berezhko, E. G., & Ellison, D. C. 1999, ApJ, 526,385 

[2] Blasi, P. 2002, Asti'opart. Phys., 16, 429 

[3] Blasi, P. 2004, Astropart. Phys., 21,45 

[4] Blasi, P, Gabici, S., & Vannoni, G. 2005, MNRAS, in press 

[5] Eichler,D. 1984, ApJ, 277,429 

[6] ElHson, D. C. 1985, J.G.R., 90, 29 

[7] ElHson, D. C, Berezhko, E. G., & Baling, M. G. 2000, ApJ, 540, 292 

[8] ElHson, D. C, DecourcheHe, A., & BaHet, J. 2004, A&A, 413, 189 

[9] Ellison, D. C. & Double, G. P 2002, Astropart. Phys., 18,213 

[10] Ellison, D. C. & Double, G. P 2004, Astropart. Phys., 22,323 

[11] Ellison, D. C, & Jones, F C, & Baring, M. G. 1999, ApJ, 512,403 

[12] Ellison, D. C, & Moebius, E., & Paschmann, G. 1990, ApJ, 352,376 

[13] Giacalone, J. and Ellison, D. C. 2000, J.G.R., 105,12,541 

[14] Jones, FC, Jokipii, J. R.& Baring, M.G. 1998, ApJ, 509,238 

[15] Jones, FC, & Ellison, D.C. 1991, Space Sci. Rev., 58,259 



